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Abstract — Magnetic Resonance Force Microscopy (MRFM) is 
an emergent technology for measuring spin-induced attonewton 
forces using a micromachined cantilever. In the interrupted Os- 
cillating Cantilever-driven Adiabatic Reversal (iOSCAR) method, 
small ensembles of electron spins are manipulated by an external 
radio frequency (RF) magnetic field to produce small periodic 
deviations in the resonant frequency of the cantilever. These 
deviations can be detected by frequency demodulation, followed 
by conventional amplitude or energy detection. In this paper, 
we develop optimal detectors for several signal models that 
have been hypothesized for measurements induced by iOSCAR 
spin manipulation. We show that two simple variants of the 
energy detector-the filtered energy detector and a hybrid filtered 
energy/amplitude/energy detector-are approximately asymptot- 
ically optimal for the Discrete-Time (D-T) random telegraph 
signal model assuming White Gaussian Noise (WGN). For the 
D-T random walk signal model, the filtered energy detector 
performs close to the optimal Likelihood Ratio Test (LRT) when 
the transition probabilities are symmetric. 



I. Introduction 

MRFM is a promising technique for dramatically improving 
the sensitivity and resolution of magnetic resonance imag- 
ing [1], [2], [3], [4]. One of the immediate goals of this 
field is to demonstrate the detection of individual electron 
spins. Successful experiments have already been performed 
that demonstrate sensitivity on the order of two electron spins 
for integration times on the order of several seconds [5]. In 
order to increase detection speed so that it is suitable for 
imaging applications, significant advances in force detection, 
spin manipulation and signal detection are required. In this 
paper, we address the topic of optimal signal detection. 

A general MRFM experiment involves the detection of 
perturbations of a thin micrometer-scale cantilever whose 
tip incorporates a submicron ferromagnet. When no electron 
spins are present, the cantilever acts as a harmonic oscillator. 
Unpaired electron spins in the sample behave like magnetic 
dipoles, exerting perturbing forces on the cantilever. Thus, 
the presence of electron spins can be detected based on 
measuring the perturbation of the cantilever position from 
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its normal oscillatory behaviour. In particular, the iOSCAR 
method uses an externally modulated RF field to manipulate 
the electron spins in such a way as to produce periodic 
forces on the oscillating cantilever [5], [6]. This results in 
small changes in the cantilever's natural frequency ujq. A laser 
interferometer measures the cantilever displacement; detection 
of these frequency shifts in the cantilever displacement signal 
identifies the presence of electron spins. 

This methodology can potentially be extended to provide 
single electron spin sensitivity. Unfortunately, there are a host 
of practical impediments to achieving this objective. Firstly, 
the spin-induced changes in c^o become extremely small at 
the single-spin level. Thus, very long integration times are 
required to detect the single spin signal. However, the integra- 
tion time is limited by spin relaxation effects that randomly 
depolarize the electron spin over time. At low temperatures, 
the relaxation effects are mitigated, which is why current 
experiments are conducted with temperatures in the millikelvin 
range. In this low temperature regime, measurements are 
sensitive to thermal noise from various sources. A major 
source of thermal noise is the heating of the cantilever by 
the laser interferometer. Spin detection methodologies must 
account for spin relaxation effects and low signal-to-noise ratio 
(SNR). 

Four signal models are presented in this paper. The first 
two are continuous-time (C-T) models while the last two are 
D-T models. There is a compelling reason for moving to D-T 
models: they are more tractable to work with in general. The 
first model is obtained by applying the classical description of 
an electron spin in a magnetic field [7]. The result is a set of 
nonlinear differential equations. The second model is derived 
using a quasi-static approximation [8]: one obtains a C-T 
random telegraph process. Neither of the two C-T models have 
finite-dimensional optimal detector implementations. In [9], 
a detector for the first model was proposed that used an 
extended Kalman Filter (KF)-like state estimator. In [10], [11], 
a hybrid Bayes/Generalized Likelihood Ratio (GLR) detector 
was developed for the C-T random telegraph model. Both have 
running times that make a real-time implementation unfeasible 
at this point. The third model is the generalized D-T equivalent 
of the C-T random telegraph, and the fourth model is a D-T 
random walk process. The optimal LRT for these last two 
models can be derived. Moreover, their running times are 
O(N) and 0(MN) respectively, where TV is the number 
of samples per observation, and the number of states in the 
D-T random walk process is 2M + 1. Surprisingly, it can 
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be shown that there exist simpler detectors, all with O(N) 
complexity, that approximate the LRT for the third model, the 
D-T random telegraph. Simulation shows that one of these 
simpler forms, the filtered energy detector, has performance 
that is comparable to the LRT for the D-T random walk under 
certain conditions. 

This paper has two main results. Firstly, the filtered en- 
ergy detector is approximately asymptotically optimal in the 
case of the symmetric D-T random telegraph model under 
the conditions of low SNR, long observation time, and the 
probability of a transition between consecutive samples (1— p) 
being small. Secondly, in the general D-T random telegraph 
model (which includes both symmetric and asymmetric tran- 
sition probabilities), a hybrid filtered energy/amplitude/energy 
detector is approximately asymptotically optimal under the 
conditions of low SNR and long observation time. The outline 
of this paper is as follows. In Section [H] we briefly review 
the iOSCAR experiment. This is followed by a discussion 
in Section HID of several iOSCAR signal models. Section [Tvl 
consists of reviewing the existing detectors that are commonly 
used, namely the amplitude and filtered energy detectors, 
and comparing them with detection schemes that we have 
developed. Simulation results are presented in Section IV1 

II. Description of the iOSCAR experiment 

A schematic description of the iOSCAR experiment at IBM 
Almaden is shown in Figure \l] In the current experiment, a 
submicron ferromagnet is placed on the tip of a cantilever 
that sits approximately 50 nanometers above a sample. In the 
presence of an applied RF field, the electron in the sample un- 
dergoes magnetic resonance if the RF field frequency matches 
the Larmor frequency. The latter is proportional to the strength 
of the cantilever tip's magnetic field. Because the tip field falls 
off rapidly with distance, only those spins that are within a thin 
resonant slice will satisfy the condition for magnetic resonance 
and interact with the cantilever. The resonant slice is located 
at a certain computable distance away from the cantilever tip. 

If the cantilever is forced into mechanical oscillation by 
positive feedback, the tip motion will cause the position of 
the resonant slice to oscillate. As the slice passes back and 
forth through an electron spin in the sample, the spin direction 
will be cyclically inverted due to an effect called adiabatic 
rapid passage [12]. The cyclic inversion is synchronous with 
the cantilever motion and affects the cantilever dynamics by 
changing the effective stiffness of the cantilever. Therefore, the 
spin-cantilever interaction can be detected by measuring small 
shifts in the period of the cantilever oscillation using a laser 
interferometer. This methodology has been successfully used 
to detect small ensembles of electron spins [5], [6]. Signal 
deconvolution of spin ensemble measurements at different 
locations above the sample and at different resonant slices 
can potentially provide single spin resolution [13]. For more 
details about iOSCAR, see [5], [6]. 

We shall briefly review the setup of the single spin- 
cantilever interaction framework proposed by Rugar et al. [7] 
and Berman et al. [8]. Consider an electron spin in a rotating 
frame that rotates at the frequency of the applied RF magnetic 
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Fig. 1. Schematic of the iOSCAR experiment. 



field Bi(t) (see Figure The effective magnetic field B e ff(t) 
in this frame is given by 



B eff {t)=B 1 (t)i + AB (t)k, 



(1) 



where i and k are the unit vectors in the x' and z directions of 
the rotating frame, B\ (t) is the amplitude of the RF magnetic 
field, Bo(t) is the amplitude of the tip magnetic field at the spin 
location, and ABo(t) = Bo(t) — wrf/7 is the off -resonance 
field amplitude. The constant 7 = 5.67T x 10 10 s _1 T _1 is the 
gyromagnetic ratio. The spins for which cjrf approximately 
equals the Larmor frequency lol = ^yBa(i) are said to be in 
magnetic resonance. Note that Bo(t) is really also a function 
of space; in our description above, we have fixed the location 
of the electron so that Bo(t) is just a function of time. Only 
electrons in a certain slice of the sample will satisfy the 
magnetic resonance condition, and the position of this slice 
is a function of the cantilever position. In the rotating frame, 
Bi(t) is a constant vector (except during the skip times which 
are dictated by the iOSCAR protocol; this will be explained in 
the next section), and AB (t) oscillates synchronously with 
the cantilever. If ABo(t) varies sufficiently slowly such that 
the adiabatic criterion 

dAB (t) 



dt 



(2) 



is satisfied, the spin can be assumed to remain aligned with 
either B e g(t) or —B e ff(t). These are the spin-lock and anti- 
spin-lock conditions, respectively. 

III. MRFM SIGNAL MODELS 

A complete analysis of the spin-cantilever interaction re- 
quires a quantum mechanical treatment. Such an analysis is 
still ongoing. Here, we shall focus on signal models that can 
be derived from a classical physics framework. A potential 
weakness of the classical approaches is that they might not 
adequately characterize the behaviour of the electron, which 
is subject to quantum effects. However, recent iOSCAR exper- 
iments have demonstrated that the key aspects of the classical 
model are valid. Experimental validation of these models will 
only be possible once successful detection of a single spin has 
been demonstrated. 

The development of the C-T random telegraph model 
(Model 2) below suggests that almost all of the information 
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Fig. 2. In the coordinate system rotating at wrf, the off-resonance field 
ABo(t), and therefore the effective field B e ff, oscillate synchronously with 
the cantilever. Under the spin-lock (anti-spin-lock) assumption, the electron 
spin aligns with (— )B e g. 



pertaining to the presence or absence of a spin is contained 
in the frequency content of the cantilever position signal z(t). 
In all of the models, we shall assume that the noise sources 
are WGN. We shall see that in the presence of a spin, z(t), 
after being frequency demodulated and translated to baseband, 
consists of an approximately periodic deterministic square 
wave and a random signal component. In the absence of the 
latter, optimal detection can be performed using a matched 
filter detector. When a random signal component is present, 
the deterministic part can be cancelled out and we are left with 
the detection of a random signal in Additive White Gaussian 
Noise (AWGN). 



accommodated by adding more second order equations similar 
to the last equation in (0, and with z%,i = 2,3, .. . used in 
the i-th additional equation in place of z. Each additional 2nd 
order equation has a different noise term F n i(t), and the z 
appearing in the first three equations of Q will be replaced 
by z + Z2 + . ■ - + z n , where n is the number of cantilever modes 
considered. Note that G ^ 0, so that when a spin is present, 
G/j, z affects the dynamics of z(t), and (0 is a nonlinear system 
of differential equations. On the other hand, when a spin is 
not present, the G/i z term vanishes, and we are left with the 
standard equation of motion for a cantilever, which is: 



Tz + kz = F n (t) 



(4) 



The observable output of the system are samples of the 
cantilever position z(t) corrupted by observation noise, which 
is assumed to be AWGN. Define t\ — iT s to be the time 
instants at which z(t) is sampled, where T s is the sampling 
interval. Model the observation noise as Wi, where Wi is 
a sequence of independent and identically distributed (i.i.d.) 
Gaussian random variables (r.v.s) with mean and variance 
a 2 . Denote the observation sample at time ti by y^. Then 
Hi = z(ti)+Wi. The detection problem for this signal model is 
as follows: given the noisy observations y = [yo, . . . , yjv-i]', 
classify the system that generated y as either: 

Ho : y generated by the no-spin system 
if i : y generated by the spin system 

In [9], we proposed a detector that uses the normal KF and 
an extended KF-like state estimator for spin detection under 
Model 1. This detector operates directly on the cantilever 
position signal. Our focus in this paper is on the detectors for 
the last two D-T models, and so we shall not make further 
mention of the dual KF detector. The interested reader is 
referred to [9] for details and results. 



A. Model 1: Classical C-T model 

The equations of the classical dynamics of a MRFM can- 
tilever interacting with a single electron spin moment are 
described in [7]. Considering only the fundamental mode 
and ignoring the positive feedback term, the interaction is 
described by: 

fix = lHy{Gz + 6 B Q ) 
fiy = i^ z Bi{t) - ~fLi x (Gz + 5B ) 
fiz = -jfiyB^t) 
mz + Tz + kz = G[i z + F n (t) (3) 

where z{t) is the position of the cantilever, z — is taken to be 
the equilibrium position, m is the cantilever's effective mass, 
and k is the cantilever spring constant. An overhead dot is 
understood to be differentiation with respect to time. The elec- 
tron spin moment is given by fl(t) = [fJ. x (t) Hy(t) l^z(t)}', 
and it is known that fj, = \fl\ = 9.28 x 1(T 24 J/T. B^t) 
is the RF signal which is known, and F n {t) is WGN which 
arises due to various noise sources in the experiment, e.g. 
background thermal noise. The above equations omit the effect 
of the higher-order modes of the cantilever. This effect can be 



B. Model 2: C-T random telegraph 

In [8], the classical C-T model (Model 1) is used to obtain 
a simpler set of equations to describe the spin-cantilever inter- 
action assuming that the cyclic adiabatic inversion condition 
(0 holds. A perturbation analysis shows that the cantilever 
position can then be described by: 

mz(t) + Tz(t) + (k + Ak)z(t) = F n (t) (5) 

Here, Ak — —^iG 2 /\Bi\. We note that the cantilever's 
natural mechanical resonance frequency is uiq = yjkjm. 
The shift in the spring constant results in a corresponding 
shift in ujo that is approximately given by — ^lu | ■ Define 
Aw* = \oJo\vG 2 /{kB 1 )\. 

With the iOSCAR protocol, B\ (t) is turned off after every 
Nskip cycles over a half-cycle duration to induce periodic 
transitions between the spin-lock and anti-spin-lock states. 
This results in Aoj alternating between the two values ±Aojq. 
By setting F n (t) = and ignoring the amplitude decay of 
z, the solution to (|3J can be approximated as a frequency- 
modulated signal: 



z(t) — Zq cos 



u t 



«(0d£ + 9 



(6) 
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where Zq is the cantilever oscillation magnitude, 9 is a random 
phase, and s(£) is a square wave that is approximately periodic 
with non-zero amplitude Awq if a spin is present and zero 
amplitude otherwise. The reason why s(£) is not periodic 
is because the oscillation period is slightly larger when the 
cantilever's natural frequency is ojq + Auiq as opposed to when 
it is luq — Aujq . However, | 1 is on the order of 10~ 6 , which 
makes s(£) approximately periodic. Thus, spin coupling (the 
presence of a spin) can be detected by frequency demodulating 
z(t) to baseband and correlating the baseband signal with a 
known square wave signal derived from B\{t). 

Unfortunately, the effects of random thermal noise and 
spin relaxation decorrelate s(£) and the square wave signal 
reference. One model for this decoherence phenomenon is sug- 
gested by the Stern-Gerlach experiment [14]: the spins main- 
tain either the spin-lock or anti-spin-lock states, but randomly 
change polarity during the course of the measurement. This 
leads to random transitions of At^o between ±Awq, which 
are assumed to have transition times distributed according to 
a Poisson process with a rate of A spin reversals/sec. Note 
that correlating the frequency demodulator output with the 
known square wave signal, as was described in the previous 
paragraph, has the effect of cancelling out the deterministic 
transitions in ujq. What remains after correlation are the 
random transitions, and as the transition times are generated 
by a Poisson process, the resultant signal takes the form of a 
so-called random telegraph process [15]. See Figure [3] 



AT. Thus, the random telegraph model is: y(t) = s(t) + w(t) 
where w(t) is AWGN with variance er^, and s(t) is a random 
telegraph signal containing only the random transitions. The 
detection problem for this model is to design a test between 
the two hypotheses: 



Hq (spin absent) 
Hi (spin present) 



y(t) = w(t) 

y(t) = s(t) + w(t) 



(7) 



for t £ [0,T]. 

A hybrid Bayes/GLR detector was previously developed 
for the C-T Random Telegraph model (Model 2) [10], [11]. 
Essentially, the detector is the LRT but with the unknown 
initial phase <fi averaged out and the Maximum Likelihood 
(ML) estimate of f and N used. The test statistic is 



mA(y) = max < In cosh 

f,N 1 



1 



'w JO 



y(t)s+(t-T,N)dt 



2*1 JO 



(s + (i;f, N)) 2 dt 



(8) 



where s + (i;r, N) is the synthesized telegraph signal having 
initial polarity = 1 (since £^[-] has been taken) and 
parametrized by f and N. As Model 2 is C-T, the parameter 
space of {r, N} is infinite-dimensional. In [10], [11], a Gibbs 
sampler was implemented to efficiently search the parameter 
space. 
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C. Model 3: Discrete-Time Random Telegraph 

Model 3 is the generalized D-T equivalent of Model 2. Here, 
we shall likewise treat {yi} 1 ^ 1 as sam pl es of the baseband 
output of the frequency demodulator and correlator. The D-T 
random telegraph signal is a D-T Markov chain, and will be 
denoted by Q, where Q e {+A 7 -A}, < i < N - 1, and Co 
is equally likely to be either ±A, A = AwJ. The transition 
probabilities of Q,i > 1 are as follows: 



P(ti\(i 



time (ms) 



V Ci = Ci-l = A 

l-p ( i = -A,( i - 1 =A 

q Ci - Ci-l = -A 

1-q Q = A, Ci_i = -A 



(9) 



Fig. 3. Top: Sample of an ideal cantilever position signal. Frequency shifts 
are not detectable by eye. Middle: Amplitude of sample RF magnetic field, 
Bi(t). It has synchronous half-cycle skips at 1 ms, 2 ms, and 3 ms for the 
creation of spin state transitions. Bottom: Ideal and noisy outputs of frequency 
demodulator under the spin presence hypothesis. It has both deterministic 
transitions due to the RF skips at 1 ms, 2 ms and 3 ms, and random ones 
due to spin relaxation. The random transitions, t, occur as a Poisson process. 
The initial polarity is <fi = 1 for this example. 



More specifically, let the baseband output of the frequency 
demodulator and correlator be denoted by y(t). Let [0, T] be 
the total measurement time period over which the correlator 
integrates the measurements, and let f — {Tj}, i — 1, . . . ,fC, 
be the time instants within this period at which random spin 
reversals occur. As r are the arrival times of a Poisson process 
with intensity A, JC is a Poisson random variable with rate 



We restrict < p, q < 1. If p = q, we say that the 
transition probabilities are symmetric, and when p ^ q, we 
shall say that they are asymmetric. For the symmetric case, 
we can match the C-T model to the D-T model by equating 
the expected number of transitions of the Poisson process to 
that of the Markov chain. This results in p = 1 — T S A. Recall 
that T s is the sampling time interval and A is the expected 
number of transitions per second. Define the signal vector ( — 
[Co, • • ■ j Ov-i]' an d me noise vector w = [wq, . . . , tojy-i]'- We 
shall model the Wi's as i.i.d. Gaussian r.v.s with mean and 
variance a 2 . The detection problem is then to decide between: 



Hq (spin absent) : y = w 
Hi (spin present) : y = £ + w 



(10) 
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D. Model 4: Brownian Motion (BM) on a sphere 



Baseband detector 



The D-T random telegraph signal model (Model 3) char- 
acterized the spin decoherence as producing random tran- 
sitions of A(jJo between ±Awq- This is a consequence of 
the assumption that the electron spin maintains either the 
spin-lock or anti-spin-lock states. An alternative model has 
been proposed which characterizes the spin as being BM on 
a sphere in R 3 [16]. Only the ^-component of the spin ft 
is measured by the cantilever: refer to the last equation of 
(0. As an approximation, we shall characterize the effect on 
Awo as producing a one-dimensional random walk confined 
to the interval I = [-AwJ, AwJ]. Discretize / into (2M + 1) 
states using a step size of s, where M £ Z and M, s > 0. 
Define Id = [—Ms, Ms] and Q to be the D-T random walk 
restricted to Id- Henceforth, we shall refer to this model 
as the D-T random walk model, keeping in mind that it 
is an approximation derived from assuming that the spin 
behaves like BM on a sphere. Associate with Q the probability 
transition matrix P, so that Pjk = P[d = (k — M)s\(,i-i = 
(j - M)s],0 < j,k < 2M, for i > 1. P is defined so that, 
at each time step, Q changes by either ±s. We shall assume 
reflecting boundary conditions in order to keep Ci in Id, and 
(o is equally likely to be either ±s. 

The detection problem is now to test (II 01 when £ is 
modelled by a random walk. Note that the D-T random walk 
model can be regarded as a multi-state generalization of the 
D-T random telegraph model. In the limit as s — > 0, M — > oo, 
the random walk converges to BM over the interval / [15]. 



IV. Frequency domain detection strategies 

The detectors considered here can be placed into three 
categories: D-T versions of existing detectors that are currently 
in use; D-T LRTs for Models 3 and 4; and approximations 
to the LRT for Model 3. The LRT is a most powerful (MP) 
test that satisfies the Neyman-Pearson criterion: it maximizes 
the probability of detection (Pd) subject to a constraint on 
the probability of false alarm (Pp) [17], which is set by 
the user. This gives us a benchmark with which to compare 
the other detectors tests for Models 3 and 4. When the 
random transition times are known, the optimal LRT is the 
matched filter, called the omniscient matched filter (MF) in this 
paper. Although unimplementable in reality, the MF detector 
provides an absolute upper bound when comparing the various 
detectors' Receiver Operating Characteristic (ROC) curves. 

The framework for the detectors is depicted in Figure |4] 
below. Note that in Figure the C-T quantity y(t) is shown 
as an input to the statistic generator; however, the detectors in 
this section operate on the sampled values {yi}. 

Define the SNR to be SNR = (lim^oo E[C,f])/a 2 , where 
d is either the D-T random telegraph or random walk process. 
I.e. the SNR is the ratio of the steady-state expected energy 
of to the noise variance. For the former process, SNR — 
A 2 /a 2 . Let the SNR in dB be SNR dB = 10 log 10 SNR. 
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Fig. 4. Baseband detector frequency demodulates the interferometric signal, 
correlates the output against a square wave p(t) whose transitions are 
synchronous with the turn-off times of the RF field B\(t), and generates 
a test statistic (e.g., accumulated squared frequency deviations), for detecting 
the presence of a spin. 



A. Amplitude, energy, filtered energy detectors 
The D-T amplitude detector is 



N 



JV-l 


Hi 




^ V 


i=a 


H 



(11) 



where rj is set to satisfy the constraint on Pp. The threshold 
j] can be empirically determined by testing il It under Ho- It 
is equivalent to the optimal LRT under the assumption that 
yi is the sum of a random constant and i.i.d. WGN. This 
assumption would be true if there were no random transitions 
in Aujo- Under such a situation, the amplitude detector is a MF 
detector. However, as the number of random transitions in yi 
increases, the performance of the amplitude detector degrades. 
An alternative test statistic is the D-T energy detector, i.e. 
the sum of the squares of the {yi} instead of the magnitude 
of the sum in (II It . As the signal and noise are assumed to 
be independent, under Hi, one would expect {yi} to have a 
higher energy on average than under Ho- This can be reliably 
detected under a sufficiently high SNR. A natural improvement 
to the energy detector is possible by pre-filtering {yi} over the 
signal bandwidth. As the signal {Ci} is baseband, a Low Pass 
Filter (LPF) is appropriate. In particular, we shall use one 
of the simplest possible LPFs: a first-order, single-pole filter 
given by 

1 - a 1 + z" 1 
Hlp(z) = — T (12) 



1 



az 



where we require \a\ < 1 for stability [18]. The time constant 
a should be chosen based on the bandwidth of the signal; 
if u c is the desired -3dB bandwidth of the filter, set a = 
(1 — siriLL> c )/ coslj c . The -3dB bandwidth depends on the 
mean number of transitions, i.e. A or 4^ in Models 2 and 
3 respectively. In practice, a bank of LPFs with different a's 
are used to perform detection. Let <Zj — yi * hi- That is, cij is 
the sequence obtained by convolving yi with hi, where hi is 
the impulse response of whatever filter we choose to use. The 
energy and filtered energy detector can then be expressed as 



N-l 



Hi 
Hn 



(13) 



6 



where for the energy detector, hi is taken to be the unit 
impulse function S[i], while for the filtered energy detector, 
hi — hup[i], the impulse response of Hup(z) in dl2l . 

We note that the running time for the amplitude, filtered 
energy, and energy detectors is O(N). 

B. Optimal LRT detectors and their approximations 

One can derive the LRTs for the D-T signal models 
(Models 3 and 4). Consider first the D-T random telegraph 
model (Model 3): define R k (S) = P(( k = S\Y k _ u . . . , Y ), 
where S G {±^1} and k > 1. Let 7V(x; p, a 2 ) 



■ ■ 



exp 



2a 2 



and 7Sl , S2 = P(Si -» S 2 ) be the 
probability that the signal Q goes from Si in the current time 
step to 52 in the next with Si, 52 G There exists a 

recursive formula for R k (S). 

R k (S) 



rVk- 



JA,S 



'Rk-iiA) 



e^ Vk -^R k _i{A) 



'fin(-i) 



+ 7 _ A ,s(l-*) 



(14) 



for k > 1 and with initial conditions i?o(^4) = ^o( — A) = 
1/2. With this, one can derive f(y; Hi), the probability density 
function (pdf) of y under i2i. Let f(y;Ho) be the pdf of y 
under Hq. Define j/ fc ) = (y k , ■ ■ ■ , yo) for k > 0. Now, 



f(y;Hi) = f(y N _i\y( N -V;Hi)- 

f(y N -2\y (N ' 3) ;Hi) ■ ■ ■ f(yi\y ; #i)/(y ; Hi} 



and 



f(Vk\v lk - l) ;H 1 ) = R k (A)M(y k ;A 7 a 2 ) 

+ R k (-A)N(y k ;-A,a 2 ), k>l 

With J15i and d!6i . the log LRT expression is: 

f(y;Hi) 



(15) 



(16) 



In A(y) = In 



W;H ) 



N-l 

^ [Rk{A)e. 

k=0 



'»+R k (-A)e 



Hi 
Ho 



(17) 



We see that, at each time step, the log LRT incorporates 
information from the present observation in exp(±Ay k / a 2 ) 
and information from the past observations in R k (±A). The 
running time of dl7> is O(N), where N is the number of 
observations. 

Under the regime of low SNR and long observation times 
(N ^> 1), the second-order expansion of < l 1 71 is approximately 
equal to the hybrid filtered energy/amplitude/energy detector: 



1 - a' 
2a 



2a 



% 2 



(18) 



where the constants C\ and Cn are given in the appendix. 
Here, a k — y k * h LP [k], i.e. the output of the observations 
convolved with the LPF in H2\ . What this means is that in the 
aforementioned regime, we expect the hybrid detector to have 



performance similar to the optimal LRT test. When p — q, the 
second-order expansion of the LRT is approximately equal to 
the filtered energy detector for values of p close to 1. See 
the appendix for more details. In light of the running times 
for the filtered energy, energy, and amplitude detectors, the 
complexity of Jl 8I > is also O(N). 

Next, consider Model 4, the D-T random walk. Define the 
vectors: 



Rk(-Ma) 



R k = R k (0) ,W k 



Rk(Ms) 



where /„(•) = Af{-;0,a 2 ). If 



fw{yk + Ms) 



fw{Vk) 



fw(Vk - Ms) 



(oi, . 



, CL ri 



and b 



(pi, . . . , b n y, define the operation a * b — (ai&i, . . . , a n b n )' 
and a - b = J^i a i^i (i- e - S • b is the dot product). Let Q = P'. 
As in^the LRT for Model 3, there exists a recursive formula 
for R k : 



Rk = 



Q(y/ k -i * fl fc -i) 

W k -i ■ R k ^i 



(19) 



for k > 1 and with the initial condition Rq — ^(cm + eM+2), 
where {e*i} are the standard basis vectors for R 2M+1 . The 
LRT for the D-T random walk (details are given in the 
appendix) can be expressed as 



N-l 



A(y) = H 



Rk-W k 



k=o e*A/+i • W k 



(20) 



The running time of the LRT for the D-T random walk is 
0(M 2 N) for a general matrix Q. If Q is tridiagonal, as is the 
case for the random walk model, the running time is 0(MN). 



V. Simulation results 

The objective in this section is to compare all of the 
detection methods discussed in this paper. The class of LRT 
detectors is optimal for their respective signal models, and 
provides a good comparison benchmark. Comparison of the 
various detectors is done using ROC curves, which is a plot 
of probability of detection (Pd) vs. probability of false alarm 
(Pf), and power curves, which is a plot of Prj vs. SNR at 
a fixed Pp. Some of the parameters used in the simulation 
of Models 3 and 4 are as follows: k = 10~ 3 N m _1 , 
luq = 2tt • 10 4 rad s"\ B x = 0.2 mT, G = 2 • 10 6 T m -1 . 
The sampling period was T s = 1 ms, and signal durations of 
T = 60 s and T = 150 s were used. The performance of 
the detectors varies as a function of T; in general, a larger T 
results in better performance. Realistic values of T are several 
orders of magnitude larger. Nevertheless, the comparative 
results obtained from using the two values of T above are 
representative of larger values. Indeed, our approximations to 
the optimal detectors improve with larger T. 
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A. D-T random telegraph model (Model 3) 

First, consider Model 3, the D-T random telegraph. Figure[5] 
depicts the simulated ROC curves at SNR = -35 dB, A = 
0.5 s _1 , and with symmetric transition probabilities (p = q). 
With T s = 1 ms, this results in p = q = 0.9995. We 
examine the matched filter, D-T random telegraph LRT (RT- 
LRT), filtered energy, hybrid, amplitude, and unfiltered energy 
detectors. The RT-LRT, filtered energy, and hybrid detector 
curves are virtually identical, which confirms our previous 
analysis. We note that the unfiltered energy and amplitude 
detectors have performance that is poorer than the RT-LRT, as 
it should be since this is the optimal detector. The unfiltered 
energy detector has the worst performance out of the five 
detector methods considered, and we shall see that this is 
almost always the case. Lastly, the omniscient MF detector 
has the best performance. Again, that is consistent with our 
expectations. We generated a power curve over a range of 
SNR under the same conditions as in Figure [6] with a fixed 
Pf = 0.1 The RT-LRT, filtered energy, and hybrid detector 
have similar performance from -30 dB to -45 dB. With this 
particular value of Pf and A, the RT-LRT, filtered energy, 
and hybrid detector perform from 5 dB to 10 dB worse than 
the MF detector. Although the amplitude detector has worse 
performance than the RT-LRT and filtered energy detector, all 
three have comparable performance at -45 dB. 






Matched filter bound 





RT LRT 


□ 


Filtered energy 




Hybrid 




Amplitude 


+ 


Unfiltered energy 



Probability of False Alarm (P ) 



Fig. 5. Simulated ROC curves for the D-T random telegraph model (Model 
3) with symmetric transition probabilities at SNR = -35 dB, T = 60 s, and 
A = 0.5 s — 1 for the omniscient matched filter, D-T random telegraph LRT, 
filtered energy, hybrid, amplitude, and unfiltered energy detectors. The RT- 
LRT is the optimal detector for this model. 



Figure Q shows the power curve generated using the bigger 
value of T = 150 s. Again, the RT-LRT, filtered energy, and 
hybrid detectors have the same performance from -30 dB 
to -45 dB. Note that the values of Prj have increased as 
compared to Figure |6] 

In the interest of space, ROC curves for a different value 
of A will not be shown. However, performance degrades as 
A increases. In any case, the curves for the RT-LRT and 
filtered energy detector are similar. Before moving on, we 
would like to present an asymmetric case where p ^ q: set 
p = 0.9998, q = 0.9992. The ROC curves are presented 





Matched filter bound 




RT LRT 


□ 


Hybrid 




Filtered energy 


■ e ■ 


Amplitude 




Unfiltered energy 



Fig. 6. Simulated power curves (Prj vs. SNR) for the D-T random telegraph 
model (Model 3) with P F fixed at 0.1 and A = 0.5 s" 1 , T = 60 s. The RT- 
LRT is the optimal detector for this model. 



- Matched tiller bound 
RT LRT 

Hybrid 

Filtered energy 
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- Unfiltered energy 



Fig. 7. Simulated power curves (Pjj vs. SNR) for the D-T random telegraph 
model (Model 3) with P F fixed at 0.1 and A = 0.5 s" 1 , T = 150 s. The 
RT-LRT is the optimal detector for this model. 



in Figure [8] There is a noticeable difference between the 
curves of the RT-LRT and filtered energy detectors. The hybrid 
detector's curve is slightly below that of the LRT, and it 
is better than that of the filtered energy detector. In fact, 
the filtered energy detector has worse performance than the 
amplitude detector. An asymmetry in p, q leads to a non-zero 
mean signal, which might be why the amplitude detector's 
performance improves. Indeed, for the D-T random telegraph 
model, lirrij—Kx, E[Q] = ^ 2 -p-q = 0-*^ f° r the values 
of p and q used here. There would therefore be larger seg- 
ments of the signal that look constant. Asymmetric transition 
probabilities can arise in some situations, e.g. experimental 
conditions or the feedback cooling of spins protocol proposed 
by Budakian [19]. 

We generated a power curve from SNR = -55 dB to -35 dB 
for the asymmetric case in Figure [5] It seems that a larger 
value of T is required when p ^ q for the hybrid filtered 
energy/amplitude/energy detector to approximate the optimal 
LRT, hence why we used T = 150 s for simulations of the 
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2M. The rest of P is: 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Probability of False Alarm (P ) 



Fig. 8. Simulated ROC curves for the D-T random telegraph model (Model 
3) with asymmetric transition probabilities (p = 0.9998, q = 0.9992) at SNR 
= -45 dB, T = 150 s. 



asymmetric random telegraph model. The hybrid detector has 
better performance than the amplitude and filtered energy 
detectors. It has performance that is comparable to the RT- 
LRT for lower SNR values. 



- Matched filter bound 

- RTLRT 
Hybrid 
Amplitude 
Filtered energy 

- Unfiltered energy 




Fig. 9. Simulated power curves (Pjj vs. SNR) for the D-T random telegraph 
model (Model 3) with P F fixed at 0.1, p = 0.9998, q = 0.9992, and T = 
150 s. The RT-LRT is the optimal detector for this model. 



p, 



K x \<i<M/2,j = i-l 

K 2 l<i<M/2,j = i+l 

0.5 M/2 < i < 3M/2, j = i - 1 or i 

Hi 3M/2 < i < 2M - l,j = i - 1 

H 2 3M/2 < i < 2M - 1, j = i + 1 



1 (21) 



In the case of M odd, the ranges for the indices i, j would 
change in an obvious way. When K\ = H 2 and K% = Hi, 
we say that the transition probabilities are symmetric, and if 
not, that they are asymmetric. In order to run the RT-LRT in 
the case of the symmetric D-T random walk, we empirically 
generate an average autocorrelation function of the random 
walk and select p (and set q — p) so that the autocorrelation 
function of the symmetric D-T random telegraph matches the 
empirical result. From this, we also obtain the optimal a for 
the LPF of the filtered energy detector. 

The ROC curves for two symmetric cases are illustrated 
in Figures [10| and El In the former, K x = K 2 = Hi = 
H 2 = 0.5, while in the latter, K x = H 2 = 0.52 and K 2 = 
Hi = 0.48. In both cases, the performance of the RW-LRT, 
RT-LRT, and filtered energy detector are all approximately 
the same, i.e. the latter two detectors are nearly optimal. 
When the transition probabilities of the D-T random walk are 
asymmetric however, as in the case of Figure El the D-T 
random walk LRT is noticeably better than the filtered energy 
detector. 
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Fig. 10. Simulated ROC curves for Model 4 at SNR = -39.9 dB, T = 60 
s, and the symmetric random walk K\ = K2 = Hi = H2 = 0.5 for the 
matched filter, D-T random walk LRT (RW-LRT), D-T random telegraph LRT 
(RT-LRT), filtered energy, amplitude, hybrid, and unfiltered energy detector. 
RW-LRT is theoretically optimal for this case. 



B. D-T random walk model (Model 4) 

For the D-T random walk model, the probability transition 
matrix P is tridiagonal. Suppose for the moment that M is 
even. Recall that the random walk Q is confined to the interval 
[—Ms, Ms]. Define the lower-quartile transition probabilities 
as Ki , K 2 and the upper-quartile transition probabilities as 
Hi, if 2- Here, we examine the performance of the detectors 
assuming the following reflecting boundary conditions: P ,i = 
1, P ,i = for i ± 1 and P 2M ,2M-i = 1, P2M,i = for i ^ 



VI. Conclusion and Discussion 

We have developed and compared optimal detectors under 
several single-spin MRFM signal models. While these models 
have to be validated through experiment, the results of this 
paper lend strong theoretical and practical support to the use 
of the simple filtered energy detector for the current MRFM 
single-spin research community. Indeed, we have shown that 
the existing baseband filtered energy detector that is in current 
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Fig. 11. Simulated ROC curves for Model 4 at SNR = -37.4 dB, T = 60 
s, and the symmetric random walk K\ = H2 = 0.52, K2 = Hi = 0.48. 
RW-LRT is theoretically optimal for this case. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Probability ot False Alarm (P p ) 



Fig. 12. Simulated ROC curves for Model 4 at SNR = -41.0 dB, T = 60 
s, and the asymmetric random walk K\ = H\ = 0.45, K2 = H2 = 0.55. 
RW-LRT is theoretically optimal for this case. 



use is approximately asymptotically optimal in the case of the 
symmetric D-T random telegraph model (Model 3) under the 
regime of low SNR, long observation times, and p close to 1 . 
The last condition can be achieved by sampling at a sufficiently 
fast rate as compared to the rate of random transitions. In the 
case of the asymmetric D-T random telegraph model, we have 
shown that a hybrid filtered energy/amplitude/energy detector 
is an approximately asymptotically optimal detector under the 
regime of low SNR and long observation times. We presented 
simulations showing that the baseband filtered energy detector 
has comparable performance with the optimal test statistic in 
the case of the symmetric D-T random walk model. We suspect 
that this similarity is not a coincidence, and work to verify this 
result analytically is ongoing. In the case of the asymmetric D- 
T random walk, the filtered energy detector does not perform 
as well as the optimal LRT We suspect that a hybrid detector 
along the lines of that formulated for the D-T random telegraph 
will perform close to the optimal. 
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Appendix I 
Derivation of the D-T random walk LRT 

Let yi be the observations, i = 0, . . . , N — 1. Put j/W = 
(yi, . . . , yo), i > and let V = {0, ±s, . . . , ±Ms} denote 
the state-space. We would like to first find f(y;H\) = 
f(y( N ~^; Hi), the density of the observations {yi} under 
Hi. From now onwards in this section, assume that the 
probabilities are conditioned on Hi unless otherwise specified. 
Define R k (S) = P(( k = S\y( k -^), S e V. Then, for k > 1, 

R k (S) = J2 P(Ck = S\( k -i = d)P(Cfc-i = d\y (k -^) 
dev 

_ Edgp P (Ck = S\( k -i = d)f w (y k _i - d)R k -i{d) 
J2iet> fwiVk-i -i)Rk-i{i) 
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which results in 



approximate solution is 

k 



Rk= (w k -i-Rk-i) 1 Q(wVi*i?fe-i) (22) fl*=/3* + ^^efc J -l/ J -,fc>0 where 



for k > 1, i?o = j(eM + ej.i+2) since (0 is equally likely to 
be either ±s. Next, for k > 1, 



/(^|y (fe - x) ) = E p (^ = d|»< fc - 1 >)/ w (s/* - d) 

dec 

=x; fl *( d )/*(»- d ) 



dec 
= Rk ■ W k 



and /(yo) = Ro 1 W^o- This leads to 

/(v (JV_1) ) = /(w-i|y (Ar - 2) ) • • ■ f(yi\yo)f(yo) 

N-l 

= ~[[Rk-W k (23) 



k=0 



The density of the observations {yi} under iJo is 

JV-l 

k=0 
N-l 

= n f^ivk) 

k=0 
N-l 

= n (^+1 • w k ) (24) 



fe=0 



and the LRT is then, using d23i and J24i : 



f(y;H ) 

N-l 

n 

fe=0 e M+i • It ft 



(25) 



Appendix II 

Approximate second-order expansion of the D-T 

RANDOM TELEGRAPH LRT AND COMPARISON TO THE 
FILTERED ENERGY TEST STATISTIC 

Let Ti (y) denote log LRT of the D-T random telegraph in 
(I17> . and T-z{y) the filtered energy detector in Jl 3i . Let us 
analyze the two test statistics under the regime of low SNR 
(|~| <§C 1) and long observation times (N 3> 1). We want 
to obtain the approximate second-order expansion of %_(y). 
Write Ti(y) ~ L\ + L 2a + L 2 b + h.o.t., where L\ are the 
1st order terms, L 2a are the 2nd order terms consisting of 
yjy k where j < k, L 2 b are the 2nd order terms of the form 
y\, and "h.o.t." are the higher order terms. Define: T k (S) = 



T k (A) 



r, for k > 0. 



R k (S)e^, S e {±A} and 9 k = T . A)+Tk{ _ A) , 
From dl4> . a recursive equation for 6 k can be derived. Its 







3=0 


Pk = 


1 - 




1 - 


r \2 


ikj = 


2(1 


— q)r k ~ j 






6fek = 


2(1 


-q) r 


1 


— r 



r k , k > 



1 - r 



1 -r 



-, < j < k - 1 
2/3 fe , k > (26) 



and r = p + q — 1. Note that p,g 6 (0, 1) =>■ |r| < 1. Define 
4y fe . Then, 



Sft 



T x ~X;{[s k (2ifc(A)-:L) + i^] 



[ Sfe (2i? fc (A)-l) + i s 2 fc ] 2 



(27) 



By solving for i?fc(^4) in terms of Ok and using J26l > in J27b . 
one can sort out the terms and obtain expressions for Li, L2 Q 



and L 2 b- Let C m = 



2-p-q 

of the mismatch in p and q. 

ii = —C m J2{l~r k )y k 



This number gives an indication 



(28) 



fc-1 



fe 3=0 



2(l- g )^._ rfc 



1 -r 



L 



2b 



A 



1-9 



y-jVk 

(29) 

(9-0(1-9) 



(1 -r) 2 



C m (2g + C m )r* + -C£r"^ 



(30) 



When p = q, C m = 0. This simplifies T\{y) considerably. 
From (EHJ-OOI, 



2 /• N-l k-1 



ft-i 



k=i 3=0 



fc=0 v ^ 



(31) 



where 7i s denotes %_ when the transition probabilities are 
symmetric. 

Next, let us obtain an expression for Tziy). For sufficiently 
large N, it can be shown that 



N-l n-l 



N-l 



Jt-j. 



k=l 3=0 



1 + a 



£»2 

k=0 



(32) 



1 — 2 

where D = n is a constant. Note that D plays no role in 

j— 1 1—1 

the performance of the test statistic. Comparing (13 11 and d32i . 
we see that they are nearly identical in form if a = 2p — 1. 
The summation of the cross-terms will be the same, but the 
coefficient of the energy term will be (1 — ^) in the case of 
T\s and (1 — -^) in the case of T%. However, the contribution 
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of —j- ^2 k y\ to T 2 is not as significant as the summation of 



the cross-terms. Now, E]J2k=o vfc H ^ ~ E [T,k=o vt H o\ 
NA 2 , and it can be shown that for large N, 



iiV-1 



JV-l fc-l 



k=l j=0 



X] XX j yjy^ H i ~ E XX!"'' J !h^ : "- 



N-l fc-l 



GA 2 (N - 1) 



where G 



q(2p-l) 



fc=l j=0 



When a = 2p-l,G 



(33) 



(2P-1) 2 



l-a(2p-l) ■ " — ^ x > w — i_( 2 p-l)2 — 

1 + i - 1. For p close to 1, G > 1, and GA 2 (iV - 1) > 



4(1— p) 1 4p 

A 2 iV. So to the first moment, the additional — ^ ^ fc y 2 to T 2 
in order to make it equal to Ti s does not represent a significant 
difference. When p « 1, we expect that the performance of 
the filtered energy detector and the LRT to be similar. 

It is possible to obtain an approximation to the second-order 
expansion of the D-T random telegraph LRT by combining 
the filtered energy, amplitude, and energy statistics. Firstly, 
for large N, the LRT is approximately 



Ti{y)=C 



(p - 



4g(l 



r)A ^ 

' < k 



Dk 



XX 

k j<k 



r k 3 y 3 Vk 



Ci 



r(l-g) 
2q{\ - r) 



X^ 



(34) 



where G = 4gj^ f-4j is constant. Let %i{y) denote the 
hybrid detector that is composed of the linear combination of 
the amplitude, filtered energy, and unfiltered energy statistics. 
We see that in the filtered energy statistic J32I . the ratio of 
the energy terms to the cross terms is yq^. For awl, this 
is roughly 1/2. The idea is to add the energy and amplitude 
statistics so that all three statistics are in the same ratio as in 
(134V So, put: 



T h (y)=T 2 (y) + 



My) 



I -a 2 
2a 



I - a' 



Gi^j/ fe + GnX 



2a 



k 



1 — a 
2a 



GnX^fc ( 35 ) 



We expect the approximation Ti(y) to have performance 
that is similar to the generalized LRT when the number of 
samples N is large. Since Th(y) is equivalent to %_{y), the 
same follows for the hybrid detector. 



